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SPECTRAL APPROACH TO D-BAR PROBLEMS 


CHRISTIAN KLEIN AND KEN MCLAUGHLIN 


Abstract. We present the first numerical approach to D-bar problems having spectral convergence for real 
analytic rapidly decreasing potentials. The proposed method starts from a formulation of the problem in 
terms of an integral equation which is solved with Fourier techniques. The singular integrand is regularized 
analytically. The resulting integral equation is approximated via a discrete system which is solved with 
Krylov methods. As an example, the D-bar problem for the Davey-Stewart son II equations is solved. The 
result is used to test direct numerical solutions of the PDE. 


1. Introduction 


It is the purpose of this paper to present an efficient numerical approach to solve D-Bar equations for a 
complex valued unknown ^ ?/), i.e., 


(i) 


di/> = ^Q(x,y)tp. 


In (1), one has (x,y) G R 2 , and we often use z = x + iy, as well as d = \ + «J0. The function Q 

is a given “potential” which we assume to be a function in the Schwartz class £>(R 2 ) of rapidly decreasing 
smooth functions in the plane. 

This equation appears in a number of different areas of analysis and applied mathematics, most promi¬ 
nently in Electrical Impedance Tomography which is of enormous importance in medical imaging, in the 
solution of completely integrable 2 + 1 dimensional partial differential equations (PDEs), as well as in two- 
dimensional orthogonal polynomials and random matrix theory (see below for details and references). 

In this paper we will primarily be interested in Complex Geometrical Optics (CGO) solutions of |l]) defined 
by the normalization at infinity, 


( 2 ) 


i[>e 


— kz 


= 1 


°'R 


where k £ C is an auxiliary parameter. 

In order to make use of the Fourier transform, we express ij) in terms of m which tends to 0 as \z\ —> oo: 
(3) m = ipe~ kz - 1 . 


This new quantity now solves 

( 4 ) dm=^e kz ~ kz Q(M+l) . 

Taking the Fourier transform defined via 

(5) F{f) = f(0 = 2. [ f{x,y)e-^ x -^dxdy , T~\f) = f(x,y) = 2. [ /(6,6)e iJlX+ ^ , 

J r2 2tt J r 2 

we find that (formally) rh solves 

( 6 ) i£m = T {Qe^ z ~ kz (m)\ + F . 

It turns out to be convenient for the numerical approach to write the equation in terms of 

(7) S(0 = £ m(£) , 
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since m is unbounded for |£| -A 0, whereas S is not. Thus the fundamental equation of study is 

(8) S{£) = -iF^Qe irz - kz Qs(0^|-^{Qe S - fez } . 

A brief gaze at Q reveals the primary obstacle for a numerical approach: the division by £, although a 
weak singularity in two dimensions, needs to be addressed. 

The approach we have developed to overcome this obstacle is to re-write the equation in a manner which 
peels off the singular behavior. The offending term appearing in (|8j) is re-written as follows: 

t ~ x G 5(o ) = T ~ l (? m) ~ am ) + ‘ f_i G g( °) ’ 

and the function G is chosen so that the following conditions are met: 

• The first term above is smooth at f = 0 (to machine precision) and vanishing (again to machine 
precision) at the boundary of the computational domain. 

• The second piece, containing all the singular terms, is evaluated via explicit analytical calculation. 
Below we choose the function G to essentially contain the ” anti-holomorphic part” of the function S , or more 
precisely a truncated Taylor series representing the anti-holomorphic part of S , multiplied by a Gaussian. 

For the numerical evaluation, the Schwartz functions are treated as essentially periodic (the computational 
domain is chosen large enough that the studied functions and their derivatives vanish to within numerical 
precision at the boundaries of the domain). The Fourier series is approximated numerically via a discrete 
Fourier series, computed via a Fast Fourier Transform (FFT). Thus the integral equation ([8]) is approximated 
via a system of linear equations for a matrix S (denoted with the same symbol in an abuse of notation) of 
finite dimension. 


In many applications, including the Davey-Stewartson equation, one often needs the solution for many 
values of the auxiliary parameter k. The approach briefly outlined above works for small values of fc|, and 
the ensuing system can be solved iteratively via GMRES [22] . 

We show, however, that for large values of k, this first approach does not work adequately, because 
the support of the desired solution is no longer distant from the boundary of the computational domain 
(for the Fourier transform variable). We develop a re-formulation of the (Fourier transformed) d equation 
based on the support properties of the solution for large k, which yields a once iterated variant of ([8]) and 
handles disparate support issues via careful shifting within the computational domain. Remarkably, this new 
equation is equivalent (up to a Fourier transform) of the iterated system of equations which Perry studies in 
([26j). and in fact can be effectively used for all values of k. 

We study the convergence of the method for concrete examples and show that the numerical approach 
has spectral convergence, i.e., that the numerical error decreases exponentially for smooth functions. The 
approach is then applied to the D-bar problem of the Davey-Stewartson (DS) II equation, and the resulting 
solution is compared to a direct numerical solution of DS II as in m- 

The paper is organized as follows: In section 2 we briefly summarize some of the main applications of 
the D-bar problem, and previous numerical approaches to its solution. In section 3 we present a short 
review on the existence theory for equation (|8|. In section 4 we discuss the numerical implementation of a 
regularization approach for |£| —> 0 of the integrand in Q. In section 5 we derive a once iterated version 
of equation ([8]) which is shown to be convenient for the construction of CGO solutions to D-bar problems, 
and use the approach of section 4 (along with careful shifting) for this equation. This approach is applied to 
the D-bar problem for the DS II system as an example in section 6. The obtained solution is compared to a 
DS II solution for the same initial data via a direct numerical integration of DS II. We add some concluding 
remarks in section 7. 

2. Applications of D-bar problems and numerical approaches 

In this section we summarize the main applications of D-bar problems, and review previous numerical 
approaches for its solution. 
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2.1. Application 1: The inverse conductivity problem. Briefly, one is given static electric measure¬ 
ments on the boundary of a domain, and the problem is to recover the conductivity in the domain. This 
is used, for example, in medical Electrical Impedance Tomography. The practical implementation involves 
tackling several different mathematical problems, one of which is precisely the d problem which is the topic 
of this paper. The mathematical formulation centers around the so-called conductivity equation, 

(9) V • ( 7 ( 2 ;, y)Vu{x, y)) = 0, (x, y) £ D . 

If one applies a voltage / on the boundary of 12, one then finds a unique solution to © together with the 
Dirichlet boundary condition 

(10) u(x, y) = f(x, y), (x, y) G <912 . 

On the other hand, if one knows the current density on the boundary (a function g(x, y)), then one can also 
find a solution to the equation ([9| together with the Neumann boundary condition 

du 

(11) i { x , y ) g ^{ x , y ) = g { x , y ), ( x , y ) edn , 

where v is the outward normal to <912. The business at hand centers on the Dirichlet-to-Neumann mapping, 
which produces an appropriate Neumann data function g given Dirichlet data /. 

The inverse problem, in a nutshell, is to start with the full knowledge of the Dirichlet-to-Neumann 
mapping, and to determine the conductivity "f(x,y) everywhere in 12. The steps, as outlined in (2,, are 
as follows: 

From the Dirichlet-to-Neumann mapping, construct the so-called scattering transform t(k). The 
steps in constructing this without direct use of the conductivity equation are outlined, for example, 
in 0 . 

Re-cast the original conductivity equation as a Schrodinger equation in the plane by writing q = 
7 _ 1 / 2 A 7 1 / 2 , and u = 7 1 / 2 u, so that 

(—A + q)u = 0 , (x,y) € 12 , 

and then extend q to the entire plane, and seek a solution u = if(x, y , k) with the following asymp¬ 
totics as \z\ —» 00 : 

e~ lkz if(x,y,k) = l + o(^j . 

The quantity if(x,y,k ) solves a D-bar equation in the variable fc, which is best expressed in terms 
of fj = e~ lkz if(x, y, k ): 

dyi = ^L 2i ( fel *- fc ^)7Z . 

47 xk 

The conductivity is extracted from g by evaluating it at k = 0: 

lim n(x, y, k) = 7 (x,y) 1/2 . 
k —>-0 

2.2. Application 2: Integrable nonlinear partial differential equations in 2 + 1 dimensions. The 

defocusing Davey-Stewartson II equation is the following nonlinear partial differential equation 


(16) 

iQt + 7 , (q X x - q yv ) = ~\q\ 2 q + m 

(17) 

Wxx + <Pyy = ^2 (|g| ) ra . 


In the above, q = q(x, y , t) = q{z , t) is complex-valued, and subscripts denote partial derivatives. Solutions 
of this equation can be interpreted (in an idealized sense) as representing the envelope of a two-dimensional 
surface wave propagating unidirectionally with a specified wave number. 

Comment on notation Here and in what follows we will write if(z, k ) even when if is not analytic in z 
or k. 


A. 

B. 

( 12 ) 

(13) 

C. 

(14) 

D. 

(15) 
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The Davey-Stewartson II equation (16 - 17) is completely integrable, possessing a Lax-pair 
(18) ip x + icr 3 ip y = I ~ 


0 q 


q 0 


(19) 


A= ( (|g| 2 ) /2 

V 


i’, 


-idq \ _ ( 0 q 

—id^d (\q\ 2 ) /2 ) U 0 


ip y + ia 3 i>. 


yy ’ 


The equation (18) is referred to, by analogy with previous 1-dimensional cases, as the spectral problem 
>o( 

( 20 ) 


associated to the DS-II equation. It can equivalently be written as 

d 0 


0 d 


^= 1 -f° « 


q 0 




The direct problem is then to seek a vector-valued solution i[) = ip(z , k ) = 
asymptotic behavior as \z\ -+ oo: 


ipi 

i>2 


with the following 


( 21 ) 

( 22 ) 


lim ipie~ Kz = 1, 

\z\ —>-oo 

lim ip 2 e~ kz = 0 . 

| z | —^OO 


As mentioned above, the quantity tj) is referred to as a CGO solution. The reflection coefficient , r(fc), is 
encoded in the sub-leading term in the asymptotic expansion of via 


(23) 




It is customary (though not necessary) to re-cast the equations (20) as a D-bar problem by defining 
(24) pi = il)ie~ kz , p 2 = ^ 2 e~ kz , 


for then p = ( ) solves the D-bar problem 

V ^2 J 


(25) 


dp = 


qe 


kz—kz 


o 1 
1 0 


l 1 


Diagonalizing the above system produces equations of the form (1). Indeed, the matrix 
eigenvalues ±1, and can be diagonalized via 


0 1 
1 0 


has 


(26) 

If one sets 
(27) 


0 1 
1 0 


1 1 
1 -1 


CT3 ' 2 


1/1 1 

1 -1 




then the entries of M = ^ ^ satisfy 


(28) 


(29) 


dMi = 


qe 


kz—kz 


-Mi, 


8M 2 = - 


qe 


kz—kz 


-Mo 


two equations of the form (jTJ) , both quantities normalized so that 


(30) 


Mj = 1 + 01-) as |z| 


oo . 
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. The reflection coefficient is obtained from (|28|)-(|29|) via 

(31) 


r 2 




- M. 


(i) 


where M-) 1 -* and m!; 1 ' 1 represent the - coefficients of M\ and M 2 as z —> oo: 


Mi = 1 + Ml 


(i) 


M 2 = 1 + 


(i) 


+ • • ' • 

z z 

As q = q{x,y,t ) evolves according to (IgJTT), the reflection coefficient evolves as well, but the evolution 
of r is quite simple: 

(32) r(k, t) = r(k, 0)e 

Better still, the recons tructio n of the potential q(x,y,t) given r(fc,0) is achieved via a D-bar problem in the 
variable k, so that in ( 33|34 ) below, d = J=. Already placing things in diagonal form, the problem is: 

(33) dMi = e kz ~ kz V{k) W l , 

(34) 8M 2 = - e kz - kz Hjc) W 2 , 


together with the normalization conditions 




M, (1) / 

f 1 \ 


(35) 

Mi 

= 1+ i + °l 

K w) 

as \k\ — t oo 



M 2 (1) 1 

(i \ 


(36) 

m 2 

= 1+ l + °( 

kW) 

as fc —> oo 


The potential q(x, y, t ) is then obtained via 

(37) q(x,y,t) = . 


2.3. Application 3: 2D orthogonal polynomials and Normal Matrix Models in Random Matrix 
Theory. The orthogonal polynomials in this example are denoted {pj = Pj(z)}^. 0 , z is the usual complex 
variable, and the orthogonality is with respect to the two-dimensional measure W(x,y)dA: 


(38) 


p k p k W(x, y)dA = S jk , j, k = 0,... 


The orthogonality condition is equivalent to the conditions contained in the statement that the following 
solid Cauchy transform decays rapidly for z —> oo: 

(39) JJ ^ W(x', y')dA' = O as z —> oo . 

Putting this together, one sees that the row vector 

(40) u~ ( Pn , HI 


solves a d-bar problem: the row vector behaves as follows for z —> oo: 


(41) 

and p, satisfies 

(42) 





0 

op = 1 1 


-W{x,y) \ 

o ) ■ 


For applications to the Normal Matrix Model in Random Matrix Theory, the sort of weight functions 
that are of interest include those of the form W = +U ( X ’A) w here the function u(x,y ) is a harmonic 

function whose growth at oo is slower than quadratic, so that the integrals are finite. While the analysis 
of orthogonal polynomials with respect to a weight supported on the real axis (or even, in some cases, 
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supported on a contour) has been carried out via Riemann-Hilbert techniques, the asymptotic analysis of 
these polynomials (as the degree grows to oo with the parameter N) has defied analysis, except in a select 
few special examples. On the other hand, their asymptotic behavior is central to a great many applications, 
from the theory of Normal Random Matrices through to approximation theory. 

2.4. Desiderata. For the first two applications described above, what is really required from the D-bar 
problem is the solution at a specific point. In the first example, the conductivity is obtained from the solution 
to the D-bar problem at k = 0, while in the second example, the reflection coefficient is obtained from the 
coefficient of = in the expansion of the solution to the D-bar problem as \z\ -A oo. (The third example, 
2D orthogonal polynomials, is more demanding, as the solution is actually required (in applications to the 
Normal Matrix Model) for all z £ C. 

In all cases, however, reformulating the D-bar equation as an integral equation shows that a starting point 
is the existence theory for a solution ip in L°°(C), which can then be promoted to a solution with greater 
regularity if desired. 

Thus, in the literature one finds existence theory in L°° , for potentials Q in relatively weak spaces, and 
one of the driving quests has been to obtain existence and uniqueness results under essentially minimal 
assumptions on the potential Q. 


2.5. Brief description of the numerical method developed by Knudsen, Mueller, and Siltanen. 

A numerical solution method for (1) was developed and implemented in |18| . The approach taken was to 
express the fundamental equation (11, normalized so that ip = l + v, with v = O (l/|z|) as \z\ -A oo, in terms 
of the operator d , which is a well-known singular integral operator in C. The integral equation takes the 
form 


(43) 


v = 1 H— 

7T 


Q{z')v 


d 2 z! 


t c 

(Compare to equation (1.3) in T8j , with k replaced by z and T(k') replaced by Q(z').) Because of the 
application to the inverse conductivity problem, they considered functions Q of compact support. This 
compact support lent itself to a periodic extension, which they showed restricts to the actual solution on the 
original support. Casting the problem in a periodic setting in turn led them to the use of the fast-Fourier 
transform in order to evaluate the discretization of the integral operator 


(44) 


T-W) = f f g{z - z')Q(z')(p(z')dz 1 dz 2 , 

J — s J — s 


where g = g(z) represents the periodic extension of the function - and the fundamental period domain is 


In discretizing and periodically extending the function 1/z, the authors handled the singularity at z = 0 
by the relatively simple rule 


(45) 


9h(z) = 


g{jh), for j £ \ 0, 


0 , 


for j = 0 . 


See [IS] for a description of the grid, and the finite integer lattice l? m . 

The authors showed that their method converges as the grid spacing tends to 0, demonstrating that the 
method is a first order method. The reason the method is first order and does not show spectral convergence 
of a Fourier approach is this simple regularization: though a Riemann integral does not change if the 
integrand is just changed in one point, the same is not true for a spectral method, in particular not for a 
discrete Fourier transform. Changing the function in one point implies that it is not continuous, and it is 
well known that the Fourier coefficients c n , n £ Z for a non-continuous function decrease as 1/n. Starting 
from the successful regularization, the authors subsequently developed a two-grid method, based on the work 
of Vainikko [23]. Note, however, that the authors consider functions of lower regularity than we, and for 
such functions spectral convergence cannot be observed even with the regularization to be discussed in the 
following sections. (Moreover, for the intended application to Electrial Impedance Tomography, the available 
data from measurements is sparse for this ill-posed problem, and so what is required is a lower order, fast 
approximation scheme.) 
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Prior to this work, there were other approaches to the numerical solution of the D-bar equation |l]) . These 
approaches are described briefly in [iSj , and we quote from that article, replacing citations with our own 
citation numbering scheme as needed (equation (1.1) in the quotation below is dv{k) = ~T(k)v(k) ): 

“The numerical solution of equation (1.1) was considered by Siltanen et al [14] for EIT in the 
numerical algorithm based on Nachman’s uniqueness proof for the 2-D inverse conductivity 
problem [23]. In [131 H511251 HB] . equation (1.1) was solved numerically by a 2-D adaptation 
of the method of product integrals [2]. The numerical solution of equation (1.1) has also been 
applied to EIT by Knudsen [19] in the numerical algorithm based on the Brown-Uhlmann 
uniqueness proof for the 2-D inverse conductivity problem [4]. A fast, direct algorithm for 
the Lippmann-Schwinger equation in two dimensions is found in [5].” 

More details can be found in Chapters 14 and 15 of the text book m- The basic approach to handling the 
inverse of the D-bar operator developed in [18] has been used in a number of different applications, see for 
example o na m ED] and references cited therein. Very recently, in [1C)I, a finite difference solver for the 
D-bar equation (which was shown to be second order) was implemented to reconstruct EIT images. 

3. Brief discussion of existence theory for equation (jT|) 

In this section we present a brief summary of known existence results for solutions to D-bar problems. 
As mentioned above, the D-bar equation |l]), re-written for the unknown m ([3]) can be re-cast as a singular 
integral equation 


(46) 


m{x,y ) = 1 + — 


Q(z')e kz '- kz 'm(z') 2 , 


an integral equation of the form 

(47) 

where 


" + 27 


(1 - ^Q.fc) TO = 1 q ±{ 1) , 


Q(z') 


r\~kz' — kz' 


(48) 


ZqM = 


2ir 


Q{z’) 


-.kz' — kz' 


Z — Z' 


f(z') 


d 2 z’ 


The existence theory for this equation has been explained in a number of places; we refer to [27] and 126] and 
references contained therein for the following description of the existence theory. (In those references, the 
authors were dedicated to establishing global existence for solutions of the Davey-Stewartson II initial value 
problem. Along the way they had to tackle the existence, uniqueness, and regularity of a coupled version of 
the above D-bar problem.) 

The operator Xq^ is first shown to be a Fredholm operator of index zero. The goal is to prove that there 
is nothing in the kernel of the operator 1 — Xq^- If this is successful, Fredholm theory implies that 1 — Xq^ 
is an invertible operator. A uniform bound on the inverse is obtained by first establishing a bound for large 
k , and then using continuity in k for the bounded subset of the k plane remaining. 

This is the original approach taken by Sung [27], who established the existence in L°°(R 2 ), with potentials 
in L 1,00 (]R 2 ) (the space of complex valued bounded and Lebesgue integrable functions on R 2 ), and then pro¬ 
ceeded to establish differentiability properties of the solution when the potential has derivatives in L 1 ’ 00 (]R 2 ). 
One particular result of his work (but not the only one) is that the reflection coefficient r{k\ 1 k 2 ) is in the 
Schwartz class Sfe li fc 2 (R 2 ) of functions if the potential Q(x,y) is in S X:V (R 2 ). 

Perry [26j iterated the integral equation one time, and thus analyzed integral operators of the form 


n kz' —kz' 


Q 


n kz"—kz" / 


■ l2 ~" d 2 z' 


He then follows the basic Fredholm theory approach outlined above, assuming now that the potential Q is 
in iJ 1 ’ 1 (R 2 ) = {/ £ L 2 (R 2 ) : V/, (1 + | • |) /(•) £ L 2 }. Note that we also use a once iterated version of the 
integral equation Q in our numerical approach, but there the iteration is done in Fourier space. 

Perry proves that if Q £ if 1,1 ) 
a unique solution in L° 


' 2 ), then the integral equation he considers (a coupled version of (46)) has 
which in addition is Holder continuous in 2 of order a for any a £ (0,1). 
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Since the main effort is in establishing that the operator (1 —Xq ^) is invertible, this implies immediately 
that our equation (46) possesses a solution in L°°(R 2 ) which is also Holder continuous with exponent a for 


oo. 


any a £ (0,1), and m actually converges to 1 as \z\ 

Moreover, if we assume that our potential is in S x>y (R 2 ), then once a solution exists which is both bounded 
and Holder continuous for z £ C, standard potential theory (see, for example, (6j pp. 53-61]) implies that the 
right hand side of (46) is Holder continuous with exponent a + 1, which clearly implies that m is infinitely 
differentiable, and possesses a complete expansion in powers of 1 for 


\Z\ -A OO. 

We summarize these facts in a Theorem, which is attributable to the above collection of references: 


Theorem 3.1. Suppose that Q £ 
expansion in powers of 4 for z 
derivatives ofm. Finally, 

(50) 


i x>y (R 2 ). Then the solution m to the d equation Q) exists, possesses an 
oo, which may be differentiated term-by-term to obtain an expansion for 


dm £ S, 


x,y\ 


2 ) 


The last conclusion, 


(50), follows from the d equation itself, and is the basis for a spectrally based 

re-written here for convenience 


numerical method. Indeed, it provides rigorous justification for equation 
(recall that S(£) = £m(£)). 


( 8 ) 


S(0 = -iJ r \Qe kz ~ kz [T-^ ( ^ 5(0 


-iT 




We close this section with a useful fact. 


Lemma 3.2. Suppose that h £ § x/! 


. Then 


TT-l 


d (h) := T-i (—T(h) 


-2 i 


is an infinitely differentiable function, possessing a complete expansion in powers of ^ as |c| 
One may see this by writing 

(51) 


-2 i. 


T~ x f -y^(£L,6)J = / h(r cos 6, r sin 9)e ir ^ x cose+vsin 9 ^~ iS d9dr , 


OO. 


which is clearly infinitely differentiable in x and y. The expansion in powers of f follows either from a 


. -=^-1 . 


stationary phase analysis of this integral, or from the alternative representation of d ( h ): 
(52) 


d \h) = 1 

7r 


z — z' 


4. Numerical algorithm for the solution of the D-bar problem 

In this section we present an algorithm for the numerical solution of the integral equation (Jsj) for small 
k. We will take k = 0, but the algorithm works for small k as well. The approach follows in principle nu, 
but uses a more sophisticated regularization of the singular integrand in (|8j) to provide a sufficiently smooth 
integrand for the numerical computation of the integral. 

4.1. Fourier approach and GMRES. It is a well known fact from Fourier analysis that the Fourier 
transform of a Schwartz function is itself in a Schwartz space. This is the analytic basis of the efficient 
implementation of so-called spectral methods, in our case discrete Fourier series, which are used in this 
paper as follows: the Schwartz functions to be studied are considered on a sufficiently large computational 
domain, 

(53) x £ L x [-n,n\, y £ L y [- 7r,7r], 

that the function and its derivatives decrease to machine precision (here roughly 10 -16 , but in practice 
limited due to rounding errors to 10 -14 or less) at the boundaries. This allows to treat Schwartz functions 
as smooth periodic functions within machine precision. In this case the Fourier coefficients decrease also 
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exponentially. Approximating the Fourier series with a discrete Fourier series (computed efficiently via an 
FFT), which corresponds in a sense to a truncation, the numerical error by neglecting the Fourier coefficients 
c n with n > N/2 will decrease exponentially with N, i.e., faster than any power of 1/N. A method being in 
this sense of infinite order is called a spectral method. We always choose the wave numbers in the FFT as 

(54) G [~N x /2 + l,N x /2]/L a , & G [~N v /2 + l,N v /2]/L v , 

where N x and N y are the number of Fourier modes in the x and y direction respectively. The values of 
N x and N y are chosen in a way that the Fourier transform of the studied Schwartz function decreases to 
machine precision as £1 and £2 approach the boundary of the computational domain. 

The discrete Fourier approach implies that the integral equation (|8j) is approximated by a finite dimensional 
system of linear equations. Writing the N x x N y matrix S as a vector with N x N y components, this system 
of equations has the form AS = b. As in m , we use GMRES [23 Krylov subspace techniques to solve 
this linear system. The advantage of GMRES is that the N x N y x N x N y matrix A does not have to be 
stored (this would make the use of parallel computers necessary for some of the values of N x , N y used here), 
just its action on a vector has to be known. The price to pay for this memory efficient implementation is 
that the equation AS 1 = b is solved iteratively with a certain tolerance that can be chosen. For the high 
precision experiments to be presented below, we typically use a tolerance of the order of machine precision 
(in practical application, this can be of course chosen larger). In practice GMRES converges rapidly in the 
studied examples. 


4.2. Regularization of the integrand. To obtain spectral convergence the goal is to provide smooth 
periodic functions for all numerical evaluations involving an FFT. Obviously this is not the case for the 
function S(f)/f in ([8]), which is not even bounded for £ —> 0, though S(£) is a smooth function, as explained 
in Theorem 3.1 To address the problem of 5(£)/£ not being in the Schwartz class, the idea is to write the 
related integral in the form 


(55) 


r 


—1 




i(S(0-G(0)+.F" 1 


^G(0 


where G(£) is chosen in a way that — G(£))/£ is (within numerical precision) a smooth rapidly decreasing 
function, and that J 7-1 (G(£)/£) can be given in explicit form analytically up to constants which can also be 
computed with spectral accuracy. If this is possible, spectral convergence can be expected since T~ l (G(£)/£) 
will be multiplied with a Schwartz function, and thus the Fourier transform on the right hand side of (55) 
can be computed with the wanted precision. 

To identify a suitable function G(£) in (55), we use the following 

Lemma 4.1. The inverse Fourier transform of exp(—1£| 2 )/£ is given by 


(56) 

In addition one has 

(57) r 
which is equivalent to 

(58) 


r- 1 1 ie-ia-1 = 


:_ (1 _ e -(x 2 +y 2 )/A 

■iy\ J 



= i(2i) r 


ni 

?n +1 




! - exp ( —— 


• ~ N n+1 


k 


£h(t) 


k =0 


k =0 


+ n+ 1) 


(zz\ "■ 

i It; ■ 


(59) 


Relation (56) can be established by direct calculation, whereas the remaining equations follow from it via 

= {-2id z ) n )~ (l-exp(-|0| 2 /4)), 


J- 


X - 


icr 




by direct calculation. 
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Remark 4.2. The term in (58) is similar to so-called ^-functions appearing in ETD schemes, see for instance 
mi and references therein. This formula is numerically problematic close to the origin (say for \z\ < 1), 


because of cancellation errors. There it is more convenient to use the Taylor expansion (58). This allows 
one to compute the function to machine precision in the whole computational domain. 


Remark 4.3. Lemma (4.1) permits the construction of a function G(£) in (55) such that (5(£) — G(£))/£ 
is a Schwartz function. A possible choice is G(£) = S(£) exp(—|£| 2 ), which implies, however, that the last 

can only be computed as a convolution, 

1 

27T 


integral on the right hand side of 
(60) T- 1 {Se~M 2 ft} = ^ I T~\S) 


z — Z‘ 


A 


1 — e 


-{z-z'){z-z')/l\ 


i2 / 

a z 


This choice of G(£) satisfies all requirements, but the computation of the convolution is not very efficient 
since the standard Fourier approach would involve functions of lower regularity than wanted (exactly the 
function exp(—|£| 2 )/£) introduced since it has an analytically known Fourier transform). Thus to obtain 
a spectral method, the convolution has to be computed directly, which is computationally expensive, since 
the integrand is an object in four real dimensions. It will be shown in this and the following section how a 
two-dimensional approach can be realized. Nonetheless we use a code implementing this approach to provide 
an independent test of other codes. 

The nature of the singularity of 5(£)/£ at £ = 0 is demonstrated by a two-dimensional Taylor expansion: 


1 5(0,0) 

r (0 - ~r 




3 


+ regular . 


The term “regular” possesses M continuous derivatives, and the first two terms clearly possess only fewer 
derivatives at £ = 0. Closer inspection shows that the first two terms can be thought of as the purely 
anti-holomorphic part of S , divided by £. Thus we choose G so as to eliminate these terms: 


(61) 


G(0 = e 


_ Mt\- 


1 <9 n 5(0) 


n—0 


n\ 


d£ r ‘ 


c 


where M is taken large enough that the Fourier coefficients of (S(£) — G (£))/£ decrease to machine precision. 
(To see this, observe that (5(£) — G(£))/£ possesses M continuous derivatives, and decays rapidly to zero for 
|£| approaching the boundary of the computational domain, and hence its Fourier coefficients will decrease 
rapidly as well). The derivatives of S in (61) are numerically computed via 


(62) 


d n S{ 0) 

<9£ n 


= I -x ) Tz n T~ x S. 


i.e., also with FFT techniques which can be done with spectral accuracy since S is a Schwartz function. 
Note, however, that for large n, the multiplication with z n in (62) will delimit the achievable numerical 
accuracy. 

To test the above approach to compute the inverse of the D-bar derivative for Schwartz functions via 
Fourier integrals, we consider the following test problem: It is straightforward to show by direct calculation 
that 


(63) 


d s 1 exp(— a(z — b){z — c)) = 




(1 — exp(— a(z — b)(z — c))) 


where a > 0, b, c are constants. In Fig. [I] we show how the difference between the numerical and the exact 
solution depends on the number of Fourier modes and the number M of terms in the Taylor expansion of S 


in (61). For simplicity we use N x = N y = N Fourier modes and L x = L y = 4. 


The left figure shows the dependence of the L°° norm of the difference between numerical and exact 
solution on the number M of terms in the Taylor series in (61), using N = 128 Fourier modes in both x 


and y directions. It can be seen in the logarithmic plot that the numerical error decreases exponentially to 
machine precision. It saturates in this example roughly for M = 11. This behavior can be understood as 
follows: as explained above, the function (S(£) — G(£))/£ possesses M continuous derivatives and decreases 
rapidly to machine precision as £ approaches the boundary of the computational domain. It is well known 
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Figure 1. Difference between the numerical solution of the D-bar inversion (63) with a = 
1/2, b = 1 and c = i. On the left we plot the dependence on the number M of terms in the 
Taylor expansion (61), with N = 128 Fourier modes in x and y. On the right we plot the 
dependence on the number of Fourier modes N in x and y 1 using M = 11 terms in the 
Taylor expansion. 


that the Fourier coefficients of such a function function decrease algebraically in the number N of Fourier 
modes as N~^ M+2 ' ) . This implies a numerical error of the same order if the series is approximated by a 
truncated series, and thus we should, and do, see an exponential decrease of the error with M. Note that 
due to the finite precision of the numerics, it does not make sense to take arbitrarily high values of M. As 


mentioned above, numerical precision is lost for large values of M in the computation of the derivatives (62) 


since the integrand is multiplied with high powers of 3 in (62) and £ in (61). Still it is possible to use values 
of M up to 20 since with larger M also the functions W n become smaller at the boundaries, and exp(—1£| 2 ) 
is small there in any case. 

In a similar way we study in the right figure of Fig. [T]the dependence of the numerical error on the number 
of Fourier modes N for fixed M = 11 (the result does not change if a higher value of M is taken). It can 
be seen that the numerical error also decreases essentially exponentially and reaches machine precision for 
iV« 2 6 . 


4.3. Algorithm for computing a solution to the D-bar equation ([8]). We present a brief summary 
for the algorithm to solve Q, with k = 0. 


(1) Introduce a discrete Fourier grid as in (53) and (54). 

(2) Compute and store the functions 


( 64 ) 


±W n = -F- 1 

n\ n\ 


r e - 


i«i 2 


using definitions (57) and (58). 


(3) Call GMRES with b = -%FQ and AS = S + iF (q J r - 1 (S'/o) ■ 

(4) To compute A(5): 

(a) Compute the relevant Taylor coefficients of S at £ = 0, using FFT techniques: 


n\ Cn = ^S( 0) = (y ) T {z n r~\S)) (0) ■ 
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(b) Form 


M 


G( ? ) = e-'f 2 , 


n —0 


(65) 


( 66 ) 


f - 1 f^G(O) = E a " W n ■ 


£ 

(c) Now compute the intermediate quantity: 

'1 


M 


n —0 




1 


M 


^ ^ CnWn 


n =0 


(d) And finally compute _4(S): 

A(S) = S + iT{Qh) . 

Note that no factorials n! that could affect the accuracy need to be computed here, since they cancel in the 


combination W n c n . Factorials just appear in the Taylor series (58). 


Remark 4.4. The computational cost of the algorithm is mainly in the FFT commands (always in two 
dimensions) in the GMRES iterations. The functions W n in (64) are computed and stored beforehand since 
they do not change during the iteration. The action of the matrix A in each GMRES iteration requires two 


FFTs. In addition the computation of the Taylor coefficients of S in (66) requires one additional FFT (not 
two, since the values are just needed for £ = 0, the second Fourier transform in (62) can be replaced by a 


summation). To sum up, per GMRES iteration 3 FFTs are necessary. Thus the computational cost is higher 
than in m (one FFT more per iteration), but the convergence of our approach is spectral instead of first 
order. When dealing with potentials of sufficient regularity (as we do), we may reach much higher precision 
with less resolution. 


5. Numerical computation of Complex Geometrical Optics solutions 


In this section we will present a numerical approach to efficiently solve the D-bar equation Q for general 
k £ C with FFT methods. Since the equation |8| defined on R 2 is approximated by a finite dimensional 
system on T 2 , the quality of the approximation will depend on the commensurability of the functions under 
consideration with the computational domain. The presence of the term e ( kz ~ kz ) in |8| implies a shift in 
Fourier space. So, for example, T( e ( kz ~ kz ) Q) will not vanish with numerical precision on the boundary of 
the computational Fourier domain for |fc| sufficiently large. Since this was the basis of the approach detailed 
in the previous section, a reformulation of ([8]) is necessary in order to apply the same techniques as before. 

5.1. Dependence on the parameter k. It is useful to re-express the fundamental equation ([8]) using an 
integral operator: 

(67) S' = /C(S) - H , 

where 


( 68 ) 


K (h) = -iT 




The operator K, is the composition of a simpler integral operator and a shift: 


(69) 


K{h)=K 0 {h)o{£ + 2ik) , 


/Co = — iT 



For sufficiently small values of fc, the numerical approach described in the previous section provides a 
spectrally accurate discretization of the equation, and the operator acting on S± can be directly inverted 

via GMRES [22]. 
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For larger values of \k\, this procedure encounters a new difficulty. Asymptotic analysis of the integral 
operator K, for \k\ large shows that it is a small norm operator, and we may invert the equation by Neumann 
series for |fc| sufficiently large. The first term of this expansion is obviously the last term on the right hand 
side of J67l): 


(70) 


S(f) « S (0) = -iJ {Qe kz ~ kz } = —iQ(£ + 2 ik) . 


Further iterations may be carried out, and so, for example, a more accu rate approximation to the solution 
is obtained by substituting S^ into the right hand side of equation (67): 


(71) 


S(0 » -i^^Qe 


D kz—kz jr-1 


+ 2 ik)^j n l — iQ(t; + 2 ik) . 


This may be simplified to the following form: 

(72) 5 (0) + 5 (1) , 
where 

(73) S w = -iT 


q 

e kz-kz-\k\ 2 /A _ g-M 2 l 

2 (z — k) 

J . 


The above asymptotic description provides a key insight into the numerical issues we are facing for k 
large. The last term is a shift of Q, and is supported in a neighborhood of —2ik, which can be very near the 
boundary of the computational domain. Moreover, the first term on the right hand side is itself the sum of 
two terms, the first also being supported in a neighborhood of —2 ik (although exponentially suppressed by 
the term e - l fe l”/ 4 ), and the second term supported near f = 0. So the solution has support in two disparate 
regimes, one near the boundary of the computational domain and the other in the middle. This behavior 
can be seen for instance in Fig. [2] 



In the original numerical approach, we are dividing by £, which is not periodic, and not very small at 
the boundary of the computational domain. This has the effect of introducing a discontinuity across the 
boundary of the computational domain, an anathema to spectral methods. 
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However, if we know beforehand that a function / has the property that its Fourier transform / has 


7^-1, 


support centered at —ik, we may compute d (/) by exploiting the shift to our advantage: 

/(£ — 2*fc)' 


(74) 


(75) 


d \f) = -2iT~ 


= -2 iT~ 


£ — 2 ik 


o (£ + 2 ik) 


= —2 ie 


kz—kz'p— 1 


/(€ - 2z_AQ 

£ — 2ifc 


Now in this last formula, the numerator, /(£ — 2 ik) is supported, by assumption, in a vicinity of £ = 0, 
and decays rapidly as |£| approaches the boundary of the computational domain. We may then apply 
the regularization approach of the previous section, but with the singularity at £ = 2 ik instead of £ = 0. 
Precisely, we take 

(76) G(0 = e-K-“l 2 £ f(0, 0 )(£ + 2ik) n , 

and then 


(77) 


d 1 (/) = -2*e 


kz—kz'p— 1 


f^-2ik)-G(0 


- 2 ie kz ~ kz T~ 


Q(0 
AC — ) 


C — 

We demonstrate the effectiveness of this shifting argument in Fig. [3] where we compute the solution to 
the example (63) for parameters leading to support of the function near the computational boundary with 
the approach of the previous section, both with and without the shifting argument It can be seen that 
the numerical solution converges as expected very slowly to the exact one without a shift since the integrand 
is not small on the boundary. A shift of the integrand in Fourier space such that its maximum is localized 
near the center of the computational domain restores rapid convergence as in the previous section. 

However, it is important to recall that the CGO solutions which we seek are actually supported (for large 
k) in two disparate regions. 



Figure 3. These are logarithmic plots of the dependence of the error in computing (631 
on N, the number of Fourier modes. The solid line shows the result of the computation 
without a shift in Fourier space, and the stars show the results with shifting. On the left, 
the parameters in (63) were chosen so that the Fourier transform of (63) is supported mid¬ 
way between the origin and the boundary of the computational domain. On the right the 
parameters were chosen so that the Fourier transform of (631 is supported very near the 
boundary of the computational domain. For these plots, the choice of k varies with TV, and 
the relevant parameter values are summarized in the Table shown in Fig. [4j 
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Fourier Modes 

L 

Left Figure k 

Right Figure k 

2 a 

4.0 

0.25(1 + 2/) 

0.5(1 +i) 

2 4 

4.0 

0.5(1+ 2/) 

1 + i 

2 b 

4.0 

(1 + 2 i) 

2 ( 1 +/) 

2 6 

4.0 

2(1 + 2 i) 

4(1+/) 

2 7 

4.0 

4(1 + 2 i) 

8 ( 1 +/) 

2 s 

4.0 

8(1 + 2 i) 

16(1+/) 


Figure 4. Parameter values for the plots shown in Fig. [3} 


5.2. Iterated integral equation. To capture the “disparate regions of support” suggested by the asymp¬ 
totic analysis, we seek a solution S of the form 


(78) 


s = m+ho(t+2ik), 


where h and / are functions supposed to vanish on the computational boundary and to have support in the 
vicinity of the origin. This leads to a pair of equations: 

(79) h = JC 0 (f)-iT(q). 

(80) / = /Co (h ° (£ + 2 ik)) o (£ + 2 ik) , 

Replacing f± in the first equation above, we have 

(81) h — /Co (ACo (h o (£ + 2 ik)) o (£ + 2 ik)) = —iT (q) . 


Note that this equation is equivalent to Perry’s integral equation (49), but that we consider it in Fourier 
space. Then we find 


(82) 


/ = /Co (h o (£ + 2ik)) o (£ + 2 ik) . 


It will turn out that in the application to the Davey Stewartson II equation, the function / is not needed, 
it just appears in the above derivation, but all the relevant information about the reflection coefficient is 
contained in the function h. 

It is useful to simplify the integral equation (81) by first noting that 
(83) 


/C 0 (h o (£ + 2 ik)) o (£ + 2 ik) = — iT I q 


I- 


£ — 2 ik 


so that the integral equation (81) becomes 
(84) h-K, of -iT ( q 


7~ 


— 2 ik 

The above simplification follows from a straightforward calculation: 


= ~i7(q) . 


(85) 

( 86 ) 

(87) 

( 88 ) 


ACo (h o + 2ik)) o (^ + 2 ik) = —iT ^<7 
= -iT | q T- 1 ( • _ ) o (^ + 2 ik) 


(; — 2 ik 


f ho (g + 2 ik) ^ 


o (^ + 2 ik) 


o (^ + 2ik) 


= -iT yi 
= = -iT I q 


kz—kzj. 7— 1 


h 


£ — 2*fc 


o (£ + 2/fc) 


CF- 


£ — 2ik 
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5.3. Algorithm to compute CGO solutions. As in section |4j we introduce a standard discretisation of 
x, y and the dual Fourier variables £ 1 , £ 2 - The role of the variable k is in a sense also dual to z as can be seen 
from equation (|8|. Thus we choose 2k on the same grid as £. The algorithm to compute the CGO solutions 
via (84) is thus as follows: 

• For given k , approximate (84) as in section [4] by discretizing 3 and £. 

• Solve the resulting linear system for the resulting vector h with GMRES. 

• The task of computing the left hand side of (|84| breaks into two parts: 


(1) Compute T 1 ^ 
compute 


h 

£-2 ik 


by the shifting and regularization argument explained in (77), and then 


(89) 


H := -iT q 




v £ — 2 ik, 

(2) Finally, compute K,${H) by using the previous regularization argument from Section 4. 

The computational cost is roughly twice the one of the algorithm of section [4] since the latter algorithm has 
to be applied twice (mainly 6 two-dimensional FFT per iteration). The advantage is that machine precision 
can be reached with much smaller values of N x , N y for all values of k with the present algorithm. 

6. Time dependence of DS II solutions 

In this section we study the time dependence of DS II solutions q(x, y, t) for given initial data q(x, y , 0). 
Parameter Convention: As noted in Subsection |5.3[ it is convenient to choose 2k on the same grid as 
£. In this Section, we make that choice, by replacing k with k/ 2, so that k and £ are dual variables and live 
on the same grid. 

6.1. Construction of solutions to the DS II equation via D-bar problems. Recall that the D-bar 
equation associated to the DS II equation is the system (28)-(29), normalized by (29). Writing S + for £M l5 
and S- for t^M^, we have two equations: 


(90) 


S ± (0 = TiJ 7 {Qe = 


T~ 


s ±{0 


=F %T 


( kz-kz I 

|Qe f I 


For each equation, we seek to represent the solution in the form 

S± = f±(0 + h ± o^ + ik) . 


It is straightforward to verify, following the same calculations from (78) to (81), that 

(91) h ± =±h , 
where h solves the equation 

(92) h — 1C 0 (/Co [h o (£ + ik)) o (£ + ik)) = —iT ( q) . 

Moreover, 

/+=/-= /Co {h o (£ + ik)) o (£ + ik) . 

For a given value of k £ C, the reflection coefficient is obtained from the computed function h via 

(93) r(k) = ih{ik) . 

This can be seen by first writing a basic representation formula for m±: 


(94) 


from which one may read off the O 


m± = — 
7 r 


< 9777,4 


>2 J 


-d z 


(95) 


,(i) 


m_i_ = — 


- term in the large z expansion: 

z 


dm± d 2 z' = (*£m±)|£_ 0 = z5±(0) . 
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Now formula (931 for the reflection coefficient follows from (951 along with (78), (91), and (82). 


To solve the DS II equation then, we compute as described above the reflection coefficient (93) and 


obtain its time dependence from (32). It can be seen that the latter has a rather simple form which 
immediately implies that |r| is constant in time which is equivalent to the existence of an infinite number 
of conserved quantities reflecting the complete integrability of the DS II equation. The most prominent 
conserved quantities are the L 2 norm of q and the energy 


(96) 


E = 


\d x q(x,y,t)\ 2 - \d v q(x,y,t)\ 2 + \q{x,y,t )\ 4 - * (</>(*, y, t) 5 


(d x 1 dy(/)(x,y,t)) 2 


It is well known that conserved quantities, the conservation of which is not implemented in the code, provide 
a test for numerical schemes. In m it was shown that the error in the numerical conservation of the DS II 
energy tends to overestimate the numerical accuracy of the solution in an L°° sense by two to three orders 
of magnitude. 

It is a remarkable fact of integrable systems that the inverse scattering techniques, here a D-bar problem, 
immediately allow to go to the final time t one is interested in. The standard approach is as follows: determine 


the reflection coefficient r{k , 0) as discussed in the previous section by solving (81) for given q(x, y, 0). Then 


solve the same equation with r(k,t) via (32) after interchanging q and r, and z and k. Thus in contrast to 
the direct numerical solution of DS II as in [131 . no intermediate time steps have to be computed. One just 


needs to solve twice integral equations of the form (81) for all considered values of k and z. 


6.2. Reconstruction of the initial data at t = 0. A first test of this approach is to recover the potential 
q(x,y,Q) from the reflection coefficient r(k, 0). To this end one fixes the parameters L x , L y in (53) and N x , 
Ny in ( |54| ). Maximal precision can only be reached for if both q(x, y, 0) and its Fourier transform decrease to 
machine precision at the boundaries of the considered domains. Thus if one wants to study the dependence 
of the difference between the computed q(x,y, 0) via the reflection coefficient and the original q(x,y, 0) on 
the parameters L x = L y = L and N x = N y = N, the pairs L, N have to be chosen such that both q(x, y, 0) 
and its Fourier transform decrease roughly to the same order of magnitude. If this is respected, then the 
numerical error decreases as expected exponentially, see the table in Fig. [5] It shows that the initial data 
can be recovered to better than 10 -13 with N = 256 Fourier modes in each direction. 


Fourier Modes 

L 

Error 

2 3 

0.7515 

7.09e-03 

2 4 

1.075 

3.1872e-04 

2 b 

1.5 

1.665e-06 

2 b 

2.1213 

1.736e-09 

2' 

3.2 

2.40e-13 

2 s 

4.2 

5.Oe-14 


Figure 5. Values of number N of Fourier modes and L for M = 11 such that both the 
reconstructed potential q(x,y,0) and the reflection coefficient decrease to roughly the same 
precision at the boundaries of the respective computational domains, and the corresponding 
errors in the reconstruction of the initial data q( x, y , 0) = exp(— x 2 — y 2 ) after solving twice 
equations of the form (81). 


6.3. DS II solutions for t > 0. It is known from the linear Schrodinger equation that Gaussian initial data 
are dispersed with time. A similar effect is to be expected in the case of the DS II solution. Thus parameters 
optimized for t = 0 will in general not provide the same resolution for t > 0. The reason is that in addition 
to the dispersive effects of DS II, there will be increasingly oscillations in the real and imaginary part of the 
DS II solution, see Fig. [G] Nonetheless, when we consider the final time t = 0.8, the computed solution still 
conserves the energy to an error of order 10~ 14 . Of course, larger values of t will require higher resolutions, 
which does not pose a principle obstacle, but might be best done on parallel computers, see the comments 
in the following section. 
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o 

'zT 


Figure 6. Solution to the DS II equation for the initial data q(x, y, 0) = exp(— x 2 — y 2 ) for 
t = 0.8; on the left the modulus of the solution, on the right the real part. 




The D-bar code can be used as a benchmark for direct solvers. For initial data in <S(R 2 ), FFT techniques 
are again the appropriate method for such solvers. The DS II equation (161 and © can be written in 
Fourier space in the form 


(97) 


i(Fq)t - (£? - £ 2 2 )^z -1 


T~ 


/c2 _ tS 
— 1 / SI S2 


CJ + S 


H\q? 


= 0 ; 


here equation (17) has been solved in Fourier space to eliminate the function p in (16). Approximating the 
Fourier transform as in the previous sections via a discrete Fourier transform, one ends up for (97) with a 
finite dimensional system of ODEs in t which are then numerically integrated with respect to t. The first 
approach along these lines appears to have been realised in [25]. More recently, in [T3], fourth order stiff 
integrators for DS II have been studied. 

The inversion of the Laplace operator in (17) has introduced the term (£ 2 — CDACi + £ 2 ) i n (97). This 
term is much better behaved than the term l/(£i +*£ 2 ) appearing in the solution of the D-bar problem, but 
it is still a problem for a spectral approach. This can be easily seen by writing ^ = |£| cos <t> and £2 = |£| sin^ 
which implies (£ 2 — £f)/(£i + Cl) = cos 2 <j> ■ Thus this factor does not have a well defined limit for £1 = £2 = 0 
and is put equal to zero there, its average value with respect to 4>. In the context of spectral accuracy as 
dicussed above, a non regular function is a problem even if there is only one non-regular point since spectral 
methods are nonlocal. 

Due to the absence of an explicitly known smooth periodic DS II solution, spectral convergence had not 
been tested in m■ The D-bar approach to DS II allows one to address this issue. Using the direct solvers 
compared in m, with data and computational parameters consistent with the example of Fig. [6j using 
Nt = 10 4 time steps, we find that the energy is conserved to the order of machine precision (the relative 
error between initial and final energy is 2e-14). The difference between this numerical solution and the one 
obtained with the D-bar code is of the order of 1CU 6 . This suggests that a regularization of the integrand in 
(97) along the lines developed in this paper for the D-bar problem will be beneficial also in the context of 
direct numerical approaches for DS II. A related approach will be discussed elsewhere. 


Remark 6.1. One can also compare the performance of the D-bar solvers in sections [4] and [5] in the context 
of DS II solutions. Using the same parameters as above for both codes, the following table shows the L°° 
norm of the difference between both solutions. Specifically, we numerically solve the inverse problem, with 
predetermined reflection coefficient ro(k±, fc 2 )e _: T( R- e ( fe )) j taking t = 0.1, in two ways: the only difference 
being in the way we computationally implement d , either using the method of Section [ 4 ] or the method 
described in this section (using the “shifting” described in Section [ 5 ]). This indicates that the more economic 
code of section [3] can be used in this context if machine precision is not needed. The reason for the relatively 
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Fourier Modes 

L 

Error 

2 b 

2.25 

7.3e-07 

2 7 

3.2 

3.4e-09 

2 s 

4.2 

5.2e-ll 


Figure 7. L°° norm of the difference between solution to DS-II at time t = 0.8 obtained 
via the D-bar approach in two ways: one way using the method of section [4] for the inversion 
of the D-bar operator, and one way using the method of section [5] for this inversion. 


high accuracy of this code is that the biggest errors are introduced there for k and 2 close to the boundaries 
of the computational domains where the reflection coefficient and the solution are smallest. Thus these errors 
contribute comparatively little to the final result. But if the CGO solution of the D-bar problem is needed 
for values of k close to the boundary with high precision, the code of section [5] has to be used. 


7. Outlook 


We have presented in this paper algorithms for the numerical computation of CGO solutions to D-bar 
problems. The presented approach is based on FFT techniques and a regularization of the functions sub¬ 
mitted to FFTs via ^-functions. It was shown that the algorithms show spectral convergence and thus allow 
computations of the solution to the D-bar problem with essentially machine precision. The computational 


cost is mainly in two-dimensional FFTs which appear in an iterative resolution of the integral equation (81) 
via GMRES. 


Whereas the convergence of the algorithm is as shown spectral, the solution of the equation (81) has to 


be computed for each considered value of k separately, which is time consuming. Thus a simple acceleration 
of the code on modern computers could be achieved by parallelizing the code in a way that the solution h 
is computed at the same time for several values of k. If higher resolutions than N x N y > 2 25 are needed, the 
corresponding FFTs are difficult to handle on a serial computer. Such resolutions can be necessary to study 
dispersive shocks , i.e., highly oscillatory regions of the DS II solution appearing in the semiclassical regime, 
see |T2] • In this case a parallelization also of the two-dimensional FFTs is necessary as in [l2j. The study of 
both types of parallelization will be the subject of further work. 

The algorithms studied in the present paper are optimized for rapidly decreasing smooth functions. For 
algebraically decreasing functions and functions with compact support being just piecewise continuous, differ¬ 
ent techniques are necessary to obtain spectral convergence. Possible approaches are multi-domain spectral 
methods as in [3] and references therein, best after introducing special coordinates as polar coordinates. This 
will be studied in an ensuing publication. 
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